######################
##  LOAD LIBRARIES  ##
######################

library(tidyverse)
library(scales)
library(modelsummary)
library(srvyr)
library(margins)


######################
##  DATA CLEANING   ##
######################

ces <- read_csv("CCES20_Common_OUTPUT.csv") %>%
  transmute(vote = case_when(CC20_410 == 1 ~ 1,
                             CC20_410 == 2 ~ 0),
            covid_self = case_when(CC20_309a_1 == 1 ~ 1,
                                   TRUE ~ 0),
            covid_knows = case_when(CC20_309a_2 == 1
                                    | CC20_309a_3 == 1 
                                    | CC20_309a_4 == 1 ~ 1,
                                    TRUE ~ 0),
            covid_death = case_when(CC20_309b_1 == 1 
                                    | CC20_309b_2 == 1
                                    | CC20_309b_3 == 1 ~ 1,
                                    TRUE ~ 0),
            age = rescale(birthyr, to = 0:1),
            gender = ifelse(gender == 2, 1, 0),
            educ = rescale(educ, to = 0:1),
            income = ifelse(faminc_new == 97, NA, faminc_new),
            income = rescale(income, to = 0:1),
            race = factor(race),
            marstat = factor(marstat),
            urbancity = factor(urbancity),
            health = rescale(CC20_309e, to = c(-1, 1)),
            income2020 = rescale(CC20_303, to = c(-1, 1)),
            economy = ifelse(CC20_302 == 6, NA, CC20_302),
            economy = rescale(economy, to = c(-1, 1)),
            pid = ifelse(pid7 == 8, NA, pid7),
            pid = rescale(-pid, to = c(-1, 1)),
            ideology = ifelse(ideo5 == 6, NA, ideo5),
            ideology = rescale(-ideology, to = c(-1, 1)),
            vote2016 = case_when(presvote16post == 1 ~ 1,
                                 presvote16post == 2 ~ -1,
                                 TRUE ~ 0),
            state = inputstate_post,
            state = case_when(state == 1 ~ "Alabama",
                              state == 2 ~ "Alaska",
                              state == 4 ~ "Arizona",
                              state == 5 ~ "Arkansas",
                              state == 6 ~ "California",
                              state == 8 ~ "Colorado",
                              state == 9 ~ "Connecticut",
                              state == 10 ~ "Delaware",
                              state == 11 ~ "District of Columbia",
                              state == 12 ~ "Florida",
                              state == 13 ~ "Georgia",
                              state == 15 ~ "Hawaii",
                              state == 16 ~ "Idaho",
                              state == 17 ~ "Illinois",
                              state == 18 ~ "Indiana",
                              state == 19 ~ "Iowa",
                              state == 20 ~ "Kansas",
                              state == 21 ~ "Kentucky",
                              state == 22 ~ "Louisiana",
                              state == 23 ~ "Maine",
                              state == 24 ~ "Maryland",
                              state == 25 ~ "Massachusetts",
                              state == 26 ~ "Michigan",
                              state == 27 ~ "Minnesota",
                              state == 28 ~ "Mississippi",
                              state == 29 ~ "Missouri",
                              state == 30 ~ "Montana",
                              state == 31 ~ "Nebraska",
                              state == 32 ~ "Nevada",
                              state == 33 ~ "New Hampshire",
                              state == 34 ~ "New Jersey",
                              state == 35 ~ "New Mexico",
                              state == 36 ~ "New York",
                              state == 37 ~ "North Carolina",
                              state == 38 ~ "North Dakota",
                              state == 39 ~ "Ohio",
                              state == 40 ~ "Oklahoma",
                              state == 41 ~ "Oregon",
                              state == 42 ~ "Pennsylvania",
                              state == 44 ~ "Rhode Island",
                              state == 45 ~ "South Carolina",
                              state == 46 ~ "South Dakota",
                              state == 47 ~ "Tennessee",
                              state == 48 ~ "Texas",
                              state == 49 ~ "Utah",
                              state == 50 ~ "Vermont",
                              state == 51 ~ "Virginia",
                              state == 53 ~ "Washington",
                              state == 54 ~ "West Virginia",
                              state == 55 ~ "Wisconsin",
                              state == 56 ~ "Wyoming"),
            weight = commonpostweight) %>% 
  drop_na()


####################
##  MAIN RESULTS  ##
####################


# Table 1 --------------------------------------------------------------------------------------------------------

models <- list(
  "(1)" = lm(vote ~ covid_self + covid_knows + covid_death,
             weights = weight,
             data = ces),
  "(2)" = lm( vote ~ covid_self + covid_knows + covid_death
              + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat),
              weights = weight,
              data = ces),
  "(3)" = lm(vote ~ covid_self + covid_knows + covid_death
             + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
             + health + income2020 + economy,
             weights = weight,
             data = ces),
  "(4)" = lm(vote ~ covid_self + covid_knows + covid_death
             + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
             + health + income2020 + economy
             + pid + ideology + vote2016,
             weights = weight,
             data = ces),
  "(5)" = lm(vote ~ covid_self + covid_knows + covid_death
             + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
             + health + income2020 + economy
             + pid + ideology + vote2016
             + factor(state),
             weights = weight,
             data = ces))

rows <- tribble(~term,                          ~"(1)",  ~"(2)",  ~"(3)",  ~"(4)",  ~"(5)",
                'Sociodemographic Controls',      'No',   'Yes',   'Yes',   'Yes',   'Yes',
                'Pandemic-related Controls',      'No',    'No',   'Yes',   'Yes',   'Yes',
                'Political Controls',             'No',    'No',    'No',   'Yes',   'Yes',
                'State Fixed-Effects',            'No',    'No',    'No',    'No',   'Yes')
attr(rows, 'position') <- 7:10

gm <- modelsummary::gof_map
gm$clean[gm$clean == 'R2 Adj.'] <- 'Adjusted R-squared'
gm$clean[gm$clean == 'Num.Obs.'] <- 'Observations'

modelsummary(models,
             vcov = "HC2",
             fmt = "%.4f",
             stars = TRUE,
             add_rows = rows,
             gof_map = gm,
             gof_omit = "AIC|BIC|Log.Lik.|F|Deviance|Sigma|Statistics|p|R2|Std. Errors",
             coef_map = c("covid_self"="Has been diagnosed with COVID-19",
                          "covid_knows"="Knows someone diagnosed with COVID-19",
                          "covid_death"="Knows someone who died from COVID-19"),
             title = "Exposure to COVID-19 and Vote Choice in the 2020 Presidential Election",
             notes = c(paste(
               "Notes: Regression estimates from linear probability models with",
               "heteroskedasticity-robust standard errors in parentheses. The dependent",
               "variable is Biden (1) versus Trump (0) in the post-election self-declared",
               "vote. Sociodemographic controls include age, gender, education level,",
               "household income, and dummies for race, marital status, and residence area.",
               "Pandemic-related controls include respondents' general health assessment,",
               "whether their household income increased or decreased in 2020, and whether",
               "they think the nation's economy got better or worse in 2020. Political",
               "controls include Democrat-Republican party identification, conservative-liberal",
               "ideology, and Trump-Clinton vote choice in 2016.")),
             output = "table1.docx")


# No-Pandemic Simulation: National-Level ---------------------------------------------------------------------------------------

means <- ces %>%
  select(covid_self, covid_knows, covid_death, 
         weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(models[["(5)"]],
           at = list(covid_self = c(0, means[[1,1]]),
                     covid_knows = c(0, means[[1,2]]),
                     covid_death = c(0, means[[1,3]])))


####################
##    APPENDIX    ##
####################


# Descriptive statistics ------------------------------------------------------------------------------------------------------

options(scipen=999)

ces_descriptives <- ces %>%
  transmute(`Vote choice (1 = Biden, 0 = Trump)` = vote,
            `Has been diagnosed with COVID-19` = covid_self,
            `Knows someone diagnosed with COVID-19` = covid_knows,
            `Knows someone who died from COVID-19` = covid_death,
            `Age (18-95)` = age,
            `Gender (female = 1)` = gender,
            `Education level` = educ,
            `Household income` = income,
            `Race: White` = ifelse(race == 1, 1, 0),
            `Race: Black or African-American` = ifelse(race == 2, 1, 0),
            `Race: Hispanic or Latino` = ifelse(race == 3, 1, 0),
            `Race: Asian or Asian-American` = ifelse(race == 4, 1, 0),
            `Race: Native American` = ifelse(race == 5, 1, 0),
            `Race: Middle Eastern` = ifelse(race == 8, 1, 0),
            `Race: Two or more races` = ifelse(race == 6, 1, 0),
            `Race: Other` = ifelse(race == 7, 1, 0),
            `Marital Status: Married` = ifelse(marstat == 1, 1, 0),
            `Marital Status: Separated` = ifelse(marstat == 2, 1, 0),
            `Marital Status: Divorced` = ifelse(marstat == 3, 1, 0),
            `Marital Status: Widowed` = ifelse(marstat == 4, 1, 0),
            `Marital Status: Never married` = ifelse(marstat == 5, 1, 0),
            `Marital Status: Domestic/civil partnership` = ifelse(marstat == 6, 1, 0),
            `Residence Area: City` = ifelse(urbancity == 1, 1, 0),
            `Residence Area: Suburb` = ifelse(urbancity == 2, 1, 0),
            `Residence Area: Town` = ifelse(urbancity == 3, 1, 0),
            `Residence Area: Rural area` = ifelse(urbancity == 4, 1, 0),
            `Residence Area: Other` = ifelse(urbancity == 5, 1, 0),
            `General health self-assessment` = health,
            `Household income increased/decreased in 2020` = income2020,
            `Nation's economy got better/worse in 2020` = economy,
            `Party identification (Republican-Democrat)` = pid,
            `Ideology (Conservative-Liberal)` = ideology,
            `2016 Vote choice (Trump-Clinton)` = vote2016)

W.Mean <- function(x) weighted.mean(x, ces$weight)
W.Median <- function(x) spatstat.geom::weighted.median(x, ces$weight)
W.SD <- function(x) radiant.data::weighted.sd(x, ces$weight)

datasummary(`Vote choice (1 = Biden, 0 = Trump)` 
            + `Has been diagnosed with COVID-19` 
            + `Knows someone diagnosed with COVID-19` 
            + `Knows someone who died from COVID-19` 
            + `Age (18-95)` 
            + `Gender (female = 1)` 
            + `Education level` 
            + `Household income` 
            + `Race: White` 
            + `Race: Black or African-American` 
            + `Race: Hispanic or Latino` 
            + `Race: Asian or Asian-American` 
            + `Race: Native American` 
            + `Race: Middle Eastern` 
            + `Race: Two or more races` 
            + `Race: Other` 
            + `Marital Status: Married`
            + `Marital Status: Separated` 
            + `Marital Status: Divorced` 
            + `Marital Status: Widowed` 
            + `Marital Status: Never married` 
            + `Marital Status: Domestic/civil partnership` 
            + `Residence Area: City`
            + `Residence Area: Suburb` 
            + `Residence Area: Town`
            + `Residence Area: Rural area` 
            + `Residence Area: Other` 
            + `General health self-assessment` 
            + `Household income increased/decreased in 2020` 
            + `Nation's economy got better/worse in 2020`
            + `Party identification (Republican-Democrat)`
            + `Ideology (Conservative-Liberal)`
            + `2016 Vote choice (Trump-Clinton)` 
            ~ (`Unique (#)` = NUnique)
            + Min * Arguments(fmt = "%.0f")
            + Max * Arguments(fmt = "%.0f")
            + (Median = W.Median)
            + (Mean = W.Mean)
            + (SD = W.SD),
            fmt = 4,
            data = ces_descriptives,
            title = "Descriptive statistics",
            output = "descriptives.docx")


# Main results using logistic regression ----------------------------------------------------------------------------

logit_models <- list(
  "(1)" = glm(vote ~ covid_self + covid_knows + covid_death,
              weights = weight, family = binomial('logit'), data = ces),
  "(2)" = glm(vote ~ covid_self + covid_knows + covid_death
              + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat),
              weights = weight, family = binomial('logit'), data = ces),
  "(3)" = glm(vote ~ covid_self + covid_knows + covid_death
              + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
              + health + income2020 + economy,
              weights = weight, family = binomial('logit'), data = ces),
  "(4)" = glm(vote ~ covid_self + covid_knows + covid_death
              + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
              + health + income2020 + economy
              + pid + ideology + vote2016,
              weights = weight, family = binomial('logit'), data = ces),
  "(5)" = glm(vote ~ covid_self + covid_knows + covid_death
              + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
              + health + income2020 + economy
              + pid + ideology + vote2016
              + factor(state), 
              weights = weight, family = binomial('logit'), data = ces))

rows <- tribble(~term,                          ~"(1)",  ~"(2)",  ~"(3)",  ~"(4)",  ~"(5)",
                'Sociodemographic Controls',      'No',   'Yes',   'Yes',   'Yes',   'Yes',
                'Pandemic-Related Controls',      'No',    'No',   'Yes',   'Yes',   'Yes',
                'Political Controls',             'No',    'No',    'No',   'Yes',   'Yes',
                'State Fixed-Effects',            'No',    'No',    'No',    'No',   'Yes')
attr(rows, 'position') <- 7:10

gm <- modelsummary::gof_map
gm$clean[gm$clean == 'Num.Obs.'] <- 'Observations'

modelsummary(logit_models,
             fmt = "%.4f",
             stars = TRUE,
             add_rows = rows,
             gof_map = gm,
             gof_omit = "Deviance Null|DF Null|DF Resid",
             coef_map = c("covid_self"="Has been diagnosed with COVID-19",
                          "covid_knows"="Knows someone diagnosed with COVID-19",
                          "covid_death"="Knows someone who died from COVID-19"),
             title = "Main results using logistic regression",
             notes = c(paste(
               "Notes: Coefficient estimates from binomial logistic regressions with standard",
               "errors in parentheses. The dependent variable is Biden (1) versus Trump (0)",
               "in the post-election self-declared vote. Sociodemographic controls include age,",
               "gender, education level, household income, and dummies for race, marital status,",
               "and residence area. Pandemic-related controls include respondents' general",
               "health assessment, whether their household income increased or decreased in",
               "2020, and whether they think the nation's economy got better or worse in 2020.",
               "Political controls include Democrat-Republican party identification,",
               "conservative-liberal ideology, and Trump-Clinton vote choice in 2016.")),
             output = "logit.docx")


# Exposure to COVID-19 and vote choice in in the four closest battlegrounds: Interactions --------------------------------------

interactions <- list(
  "State heterogeneous effects" = lm(vote ~ covid_self + covid_knows + covid_death
                                     + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                                     + health + income2020 + economy
                                     + pid + ideology + vote2016
                                     + battleground 
                                     + covid_self * battleground
                                     + covid_knows * battleground
                                     + covid_death * battleground, 
                                     weights = weight, 
                                     data = ces %>%
                                       mutate(battleground = case_when(state == "Arizona" ~ "Arizona",
                                                                       state == "Pennsylvania" ~ "Pennsylvania",
                                                                       state == "Wisconsin" ~ "Wisconsin",
                                                                       state == "Georgia" ~ "Georgia",
                                                                       TRUE ~ "0"))))

rows <- tribble(~term,  ~"State heterogeneous effects",
                'Sociodemographic Controls',     'Yes',
                'Pandemic-Related Controls',     'Yes',
                'Political Controls',            'Yes',
                'State Fixed-Effects',           'No')
attr(rows, 'position') <- 20:23

gm <- modelsummary::gof_map
gm$clean[gm$clean == 'R2 Adj.'] <- 'Adjusted R-squared'
gm$clean[gm$clean == 'Num.Obs.'] <- 'Observations'

cm <- c("covid_self"="Has been diagnosed with COVID-19",
        "covid_knows"="Knows someone diagnosed with COVID-19",
        "covid_death"="Knows someone who died from COVID-19",
        "battlegroundGeorgia"="Georgia",
        "covid_self:battlegroundGeorgia"="Has been diagnosed with COVID-19 × Georgia",
        "covid_knows:battlegroundGeorgia"="Knows someone diagnosed with COVID-19 × Georgia",
        "covid_death:battlegroundGeorgia"="Knows someone who died from COVID-19 × Georgia",
        "battlegroundArizona"="Arizona",
        "covid_self:battlegroundArizona"="Has been diagnosed with COVID-19 × Arizona",
        "covid_knows:battlegroundArizona"="Knows someone diagnosed with COVID-19 × Arizona",
        "covid_death:battlegroundArizona"="Knows someone who died from COVID-19 × Arizona",
        "battlegroundWisconsin"="Wisconsin",
        "covid_self:battlegroundWisconsin"="Has been diagnosed with COVID-19 × Wisconsin",
        "covid_knows:battlegroundWisconsin"="Knows someone diagnosed with COVID-19 × Wisconsin",
        "covid_death:battlegroundWisconsin"="Knows someone who died from COVID-19 × Wisconsin",
        "battlegroundPennsylvania"="Pennsylvania",
        "covid_self:battlegroundPennsylvania"="Has been diagnosed with COVID-19 × Pennsylvania",
        "covid_knows:battlegroundPennsylvania"="Knows someone diagnosed with COVID-19 × Pennsylvania",
        "covid_death:battlegroundPennsylvania"="Knows someone who died from COVID-19 × Pennsylvania")

modelsummary(interactions,
             vcov = "HC2",
             fmt = "%.4f",
             stars = TRUE,
             add_rows = rows,
             gof_map = gm,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             gof_omit = "AIC|BIC|Log.Lik.|F|Deviance|Sigma|Statistics|p|R2|Std. Errors",
             coef_map = cm,
             title = "Exposure to COVID-19 and vote choice in in the four closest battlegrounds: Interactions ",
             notes = c(paste(
               "Notes: Regression estimates from linear probability models with",
               "heteroskedasticity-robust standard errors in parentheses.",
               "The dependent variable is Biden (1) versus Trump (0) in the",
               "post-election self-declared vote. Sociodemographic controls include",
               "age, gender, education level, household income, and dummies for race,",
               "marital status, and residence area. Pandemic-related controls include",
               "respondents' general health assessment, whether their household income",
               "increased or decreased in 2020, and whether they think the nation's",
               "economy got better or worse in 2020. Political controls include",
               "Democrat-Republican party identification, conservative-liberal",
               "ideology, and Trump-Clinton vote choice in 2016.")),
             output = "interactions.docx")


# Exposure to COVID-19 and vote choice in the four closest battlegrounds: Aggregated states --------------------------------------

battlegrounds <- list(
  "Four Closest States" = lm(vote ~ covid_self + covid_knows + covid_death 
                             + health + income2020 + economy
                             + factor(race) + educ + age + gender 
                             + income + factor(urbancity) + factor(marstat)
                             + pid + ideology + vote2016
                             + factor(state), 
                             weights = weight, 
                             data = subset(ces, 
                                           state %in% c("Pennsylvania", "Wisconsin", 
                                                        "Arizona", "Georgia"))))

gm <- modelsummary::gof_map
gm$clean[gm$clean == 'R2 Adj.'] <- 'Adjusted R-squared'
gm$clean[gm$clean == 'Num.Obs.'] <- 'Observations'

rows <- tribble(~term,                          ~"Four Battlegrounds",
                'Sociodemographic Controls',     'Yes', 
                'Pandemic-Related Controls',      'Yes',
                'Political Controls',           'Yes',
                'State Fixed-Effects',           'Yes')
attr(rows, 'position') <- 7:10

modelsummary(battlegrounds,
             vcov = "HC2",
             fmt = "%.3f",
             stars = TRUE,
             add_rows = rows,
             gof_map = gm,
             gof_omit = "AIC|BIC|Log.Lik.|F|Deviance|Sigma|Statistics|p|R2|Std. Errors",
             coef_map = c("covid_self"="Has been diagnosed with COVID-19",
                          "covid_knows"="Knows someone diagnosed with COVID-19",
                          "covid_death"="Knows someone who died from COVID-19"),
             title = "Exposure to COVID-19 and vote choice in the four closest battlegrounds: Aggregated states",
             notes = c(paste(
               "Notes: Regression estimates from linear probability models with",
               "heteroskedasticity-robust standard errors in parentheses. The dependent",
               "variable is Biden (1) versus Trump (0) in the post-election self-declared",
               "vote. Sociodemographic controls include age, gender, education level,",
               "household income, and dummies for race, marital status, and residence area.",
               "Pandemic-related controls include respondents' general health assessment,",
               "whether their household income increased or decreased in 2020, and whether",
               "they think the nation's economy got better or worse in 2020. Political",
               "controls include Democrat-Republican party identification, conservative-liberal",
               "ideology, and Trump-Clinton vote choice in 2016.")),
             output = "four_battlegrounds.docx")


# Exposure to COVID-19 and vote choice in the four closest battlegrounds: Individual states

individual_states <- list(
  "Georgia" = lm(vote ~ covid_self + covid_knows + covid_death
                 + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                 + health + income2020 + economy
                 + pid + ideology + vote2016,
                 weights = weight, 
                 data = subset(ces, state == "Georgia")),
  "Arizona" = lm(vote ~ covid_self + covid_knows + covid_death
                 + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                 + health + income2020 + economy
                 + pid + ideology + vote2016,
                 weights = weight,
                 data = subset(ces, state == "Arizona")),
  "Wisconsin" = lm(vote ~ covid_self + covid_knows + covid_death
                   + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                   + health + income2020 + economy
                   + pid + ideology + vote2016,
                   weights = weight,
                   data = subset(ces, state == "Wisconsin")),
  "Pennsylvania" = lm(vote ~ covid_self + covid_knows + covid_death
                      + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                      + health + income2020 + economy
                      + pid + ideology + vote2016,
                      weights = weight,
                      data = subset(ces, state == "Pennsylvania")))

rows <- tribble(~term,                     ~"Georgia",  ~"Arizona",  ~"Wisconsin", ~"Pennsylvania",
                'Sociodemographic Controls',    'Yes',       'Yes',         'Yes',           'Yes', 
                'Pandemic-Related Controls',    'Yes',       'Yes',         'Yes',           'Yes',
                'Political Controls',           'Yes',       'Yes',         'Yes',           'Yes')
attr(rows, 'position') <- 7:10

gm <- modelsummary::gof_map
gm$clean[gm$clean == 'R2 Adj.'] <- 'Adjusted R-squared'
gm$clean[gm$clean == 'Num.Obs.'] <- 'Observations'

modelsummary(individual_states,
             vcov = "HC1",
             fmt = "%.2f",
             stars = TRUE,
             add_rows = rows,
             gof_map = gm,
             gof_omit = "AIC|BIC|Log.Lik.|F|Deviance|Sigma|Statistics|p|R2|Std. Errors",
             coef_map = c("covid_self"="Has been diagnosed with COVID-19",
                          "covid_knows"="Knows someone diagnosed with COVID-19",
                          "covid_death"="Knows someone who died from COVID-19"),
             title = "Exposure to COVID-19 and vote choice in the four closest battlegrounds: Individual states",
             notes = c(paste(
               "Notes: Regression estimates from linear probability models with",
               "heteroskedasticity-robust standard errors in parentheses. The dependent",
               "variable is Biden (1) versus Trump (0) in the post-election self-declared",
               "vote. Sociodemographic controls include age, gender, education level,",
               "household income, and dummies for race, marital status, and residence area.",
               "Pandemic-related controls include respondents' general health assessment,",
               "whether their household income increased or decreased in 2020, and whether",
               "they think the nation's economy got better or worse in 2020. Political",
               "controls include Democrat-Republican party identification, conservative-liberal",
               "ideology, and Trump-Clinton vote choice in 2016.")),
             output = "individual_states.docx")


# Exposure to COVID-19 and vote choice among the overall electorate, Democrats, Independents, and Republicans ------------------------------

pid_models <- list(
  "Overall (1)" = lm(vote ~ covid_self + covid_knows + covid_death
                     + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                     + health + income2020 + economy
                     + pid + ideology + vote2016
                     + factor(state),
                     weights = weight,
                     data = ces),
  "Democrats (2)" = lm(vote ~ covid_self + covid_knows + covid_death
                       + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                       + health + income2020 + economy
                       + pid + ideology + vote2016
                       + factor(state),
                       weights = weight,
                       data = subset(ces, pid > .5)),
  "Independents (3)" = lm(vote ~ covid_self + covid_knows + covid_death
                          + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                          + health + income2020 + economy
                          + pid + ideology + vote2016
                          + factor(state),
                          weights = weight,
                          data = subset(ces, between(pid, -.5, .5))),
  "Republicans (4)" = lm(vote ~ covid_self + covid_knows + covid_death
                         + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                         + health + income2020 + economy
                         + pid + ideology + vote2016
                         + factor(state),
                         weights = weight,
                         data = subset(ces, pid < -.5)))

rows <- tribble(~term,                           ~"(1)", ~"(2)", ~"(3)", ~"(4)",
                'Sociodemographic Controls',      'Yes',  'Yes',  'Yes',  'Yes',
                'Pandemic-Related Controls',      'Yes',  'Yes',  'Yes',  'Yes',
                'Political Controls',             'Yes',  'Yes',  'Yes',  'Yes',
                'State Fixed-Effects',            'Yes',  'Yes',  'Yes',  'Yes')
attr(rows, 'position') <- 7:10

gm <- modelsummary::gof_map
gm$clean[gm$clean == 'R2 Adj.'] <- 'Adjusted R-squared'
gm$clean[gm$clean == 'Num.Obs.'] <- 'Observations'

modelsummary(pid_models,
             vcov = "HC2",
             fmt = 3,
             stars = TRUE,
             add_rows = rows,
             gof_map = gm,
             gof_omit = "AIC|BIC|Log.Lik.|F|Deviance|Sigma|Statistics|p|R2|Std. Errors",
             coef_map = c("covid_self"="Has been diagnosed with COVID-19",
                          "covid_knows"="Knows someone diagnosed with COVID-19",
                          "covid_death"="Knows someone who died from COVID-19"),
             title = "Exposure to COVID-19 and vote choice among the overall electorate, Democrats, Independents, and Republicans",
             notes = c(paste(
               "Notes: Regression estimates from linear probability models with",
               "heteroskedasticity-robust standard errors in parentheses. The dependent",
               "variable is Biden (1) versus Trump (0) in the post-election self-declared",
               "vote. Sociodemographic controls include age, gender, education level,",
               "household income, and dummies for race, marital status, and residence area.",
               "Pandemic-related controls include respondents' general health assessment,",
               "whether their household income increased or decreased in 2020, and whether",
               "they think the nation's economy got better or worse in 2020. Political",
               "controls include Democrat-Republican party identification, conservative-liberal",
               "ideology, and Trump-Clinton vote choice in 2016. The first column corresponds",
               "to the fully specified model (Column 5 of Table 1).")),
             output = "pid_results.docx")


# Exposure to COVID-19 and vote choice in North Carolina, Nevada, Michigan, and Florida -------------------------------------------------

individual_states_other <- list(
  "North Carolina" = lm(vote ~ covid_self + covid_knows + covid_death
                        + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                        + health + income2020 + economy
                        + pid + ideology + vote2016,
                        weights = weight,
                        data = subset(ces, state == "North Carolina")),
  "Nevada" = lm(vote ~ covid_self + covid_knows + covid_death
                + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                + health + income2020 + economy
                + pid + ideology + vote2016,
                weights = weight,
                data = subset(ces, state == "Nevada")),
  "Michigan" = lm(vote ~ covid_self + covid_knows + covid_death
                  + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                  + health + income2020 + economy
                  + pid + ideology + vote2016,
                  weights = weight,
                  data = subset(ces, state == "Michigan")),
  "Florida" = lm(vote ~ covid_self + covid_knows + covid_death
                 + factor(race) + educ + age + gender + income + factor(urbancity) + factor(marstat)
                 + health + income2020 + economy
                 + pid + ideology + vote2016,
                 weights = weight,
                 data = subset(ces, state == "Florida")))

rows <- tribble(~term,                 ~"North Carolina",  ~"Nevada",  ~"Michigan",  ~"Florida", 
                'Sociodemographic Controls',       'Yes',      'Yes',        'Yes',       'Yes',  
                'Pandemic-Related Controls',       'Yes',      'Yes',        'Yes',       'Yes', 
                'Political Controls',              'Yes',      'Yes',        'Yes',       'Yes')
attr(rows, 'position') <- 7:10

gm <- modelsummary::gof_map
gm$clean[gm$clean == 'R2 Adj.'] <- 'Adjusted R-squared'
gm$clean[gm$clean == 'Num.Obs.'] <- 'Observations'

modelsummary(individual_states_other,
             fmt = 2,
             vcov = "HC1",
             stars = TRUE,
             add_rows = rows,
             gof_map = gm,
             gof_omit = "AIC|BIC|Log.Lik.|F|Deviance|Sigma|Statistics|p|R2|Std. Errors",
             coef_map = c("covid_self"="Has been diagnosed with COVID-19",
                          "covid_knows"="Knows someone diagnosed with COVID-19",
                          "covid_death"="Knows someone who died from COVID-19"),
             title = "Exposure to COVID-19 and vote choice in North Carolina, Nevada, Michigan, and Florida",
             notes = c(paste(
               "Notes: Regression estimates from linear probability models with",
               "heteroskedasticity-robust standard errors in parentheses. The dependent",
               "variable is Biden (1) versus Trump (0) in the post-election self-declared",
               "vote. Sociodemographic controls include age, gender, education level,",
               "household income, and dummies for race, marital status, and residence area.",
               "Pandemic-related controls include respondents' general health assessment,",
               "whether their household income increased or decreased in 2020, and whether",
               "they think the nation's economy got better or worse in 2020. Political",
               "controls include Democrat-Republican party identification, conservative-liberal",
               "ideology, and Trump-Clinton vote choice in 2016.")),
             output = "individual_states_other.docx")

# No-Pandemic Simulation: State-Level -------------------------------------------------------------------------------

four_closest_states <- ces %>%
  filter(state %in% c("Pennsylvania", "Wisconsin", 
                      "Arizona", "Georgia")) %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(battlegrounds[["Four Closest States"]],
           at = list(covid_self = c(0, four_closest_states[1,1]),
                     covid_knows = c(0, four_closest_states[1,2]),
                     covid_death = c(0, four_closest_states[1,3])))

Pennsylvania <- ces %>%
  filter(state == "Pennsylvania") %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(individual_states[["Pennsylvania"]], 
           at = list(covid_self = c(0, Pennsylvania[1,1]),
                     covid_knows = c(0, Pennsylvania[1,2]),
                     covid_death = c(0, Pennsylvania[1,3])))

Georgia <- ces %>%
  filter(state == "Georgia") %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(individual_states[["Georgia"]], 
           at = list(covid_self = c(0, Georgia[1,1]),
                     covid_knows = c(0, Georgia[1,2]),
                     covid_death = c(0, Georgia[1,3])))

Arizona <- ces %>%
  filter(state == "Arizona") %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(individual_states[["Arizona"]], 
           at = list(covid_self = c(0, Arizona[1,1]),
                     covid_knows = c(0, Arizona[1,2]),
                     covid_death = c(0, Arizona[1,3])))

Wisconsin <- ces %>%
  filter(state == "Wisconsin") %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(individual_states[["Wisconsin"]], 
           at = list(covid_self = c(0, Wisconsin[1,1]),
                     covid_knows = c(0, Wisconsin[1,2]),
                     covid_death = c(0, Wisconsin[1,3])))

Nevada <- ces %>%
  filter(state == "Nevada") %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(individual_states_other[["Nevada"]], 
           at = list(covid_self = c(0, Nevada[1,1]),
                     covid_knows = c(0, Nevada[1,2]),
                     covid_death = c(0, Nevada[1,3])))

Florida <- ces %>%
  filter(state == "Florida") %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(individual_states_other[["Florida"]], 
           at = list(covid_self = c(0, Florida[1,1]),
                     covid_knows = c(0, Florida[1,2]),
                     covid_death = c(0, Florida[1,3])))

NC <- ces %>%
  filter(state == "North Carolina") %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(individual_states_other[["North Carolina"]], 
           at = list(covid_self = c(0, NC[1,1]),
                     covid_knows = c(0, NC[1,2]),
                     covid_death = c(0, NC[1,3])))

Michigan <- ces %>%
  filter(state == "Michigan") %>%
  select(covid_self, covid_knows, covid_death, weight) %>%
  as_survey(weights = c(weight)) %>%
  summarize(across(1:3, survey_mean, vartype = NULL))

prediction(individual_states_other[["Michigan"]], 
           at = list(covid_self = c(0, Michigan[1,1]),
                     covid_knows = c(0, Michigan[1,2]),
                     covid_death = c(0, Michigan[1,3])))


